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Abstract. The evolution of an initially prepared distribution of micron sized colloidal particles, trapped at 
a fluid interface and under the action of their mutual capillary attraction, is analyzed by using Brownian 
dynamics simulations. At a separation A given by the capillary length of typically 1 mm, the distance 
dependence of this attraction exhibits a crossover from a logarithmic decay, formally analogous to two- 
dimensional gravity, to an exponential decay. We discuss in detail the adaption of a particle-mesh algorithm, 
as used in cosmological simulations to study structure formation due to gravitational collapse, to the present 
colloidal problem. These simulations confirm the predictions, as far as available, of a mean-field theory 
developed previously for this problem. The evolution is monitored by quantitative characteristics which are 
particularly sensitive to the formation of highly inhomogeneous structures. Upon increasing A the dynamics 
show a smooth transition from the spinodal decomposition expected for a simple fluid with short-ranged 
attraction to the self-gravitational collapse scenario. 



PACS. 82.70.Dd - 47.11. Mn - 05.40. Jc 

1 Introduction 

The investigation of systems in low dimensions offers a rich 
spectrum of phenomena which are rather distinct from 
three-dimensional (3D) bulk properties. This is the case 
because structural properties of condensed matter depend 
' sensitively on the spatial dimension. For example, in two- 
dimensional (2D) systems the solid-liquid phase transi- 
tion is replaced by a solid-hexatic-liquid transition [TJ[2J 
' [3] . Micron sized colloidal particles with radii i?o trapped 
at a fluid interface form a monolayer and constitute an 
• easily accessible mesoscopic model system which exhibits 
, many characteristics of truly 2D systems. This system can 
■ actually have fluid and crystalline phases [3], as well as a 
hexatic phase [5] . Applications may involve controlled self 
, assembly [6] and structure formation [ZEE], exploiting 
' the specific features of the interactions between colloidal 
particles. 

Specifically, the focus of the present study is 2D parti- 
cle aggregation due to interactions induced by interfacial 
deformation, i.e., capillary attraction (see, e.g., Refs. [TU1 
[TT] and references therein) . We consider a collection of col- 
loidal particles at the interface between two fluid phases, 
each colloid being acted upon vertically by an external 
force / (due to, e.g., gravity, an external electric field, an 
optical tweezer, etc.; see Fig.Q]and Ref. }12)b 

In the limit of small interfacial deformations or, al- 
ternatively, if the interparticle separation d between the 
centers of the colloids is large enough, the effective capil- 




f - lower fluid phase 



Fig. 1. a) A single colloid adsorbed at a fluid interface (dot- 
ted line) with no external forces present, b) Identical particles 
trapped at the interface between two fluids and deforming it 
due to an external force / acting on each of them [T2] Upwards 
pointing forces are taken to be positive. In all cases the inter- 
face meets the surfaces of the colloids under a fixed contact 
angle (9. 

lary attraction is described by the pair potential [TBI 141 

[JMTHEEMIS] 

WO = -V K (~) , V :=^; (1) 



A := sj^j / {gAp) is the capillary length, which determines 
the range of the interaction and is given in terms of the 
gravitational acceleration g, the surface tension 7 of the 
interface, and the mass density difference Ap between the 
two fluid phases. The prefactor Vq sets the strength of the 
interaction and Kq denotes a modified Bcsscl function. 
Due to V cap (d <C A) ~ hid, in that range the potential 
is non-integrable (in the sense of equilibrium statistical 
mechanics, i.e., leading to a hyperextensive scaling of the 



2 



J. Bleibel, A. Dommguez, M. Oettel, S. Dietrich: Collective dynamics of colloids at fluid interfaces 



energy), while for A finite V ca , p (d A) ~ cxp(— d/X). This 
functional dependence is analogous to screened electro- 
statics or gravity in 2D, with A playing the role of the 
screening length and / the role of the electric charge or 
the gravitational mass. The difference with electrostatics 
is that here equal (different) charges attract (repel) each 
other, and the difference with gravity is that here one can 
have negative as well as positive masses. Typically, A is 
of the order of a few millimeters so that A 3> Rq. This 
implies that a single particle interacts with many neigh- 
bors through a non-integrable potential and, in this sense 
and range, the capillary attraction formally qualifies as 
a long-ranged interaction. Long-ranged interactions en- 
joy significant current interest (see, e.g., Ref. [3U] and ref- 
erences therein) in many fields, such as self-gravitating 
fluids in cosmology, two-dimensional vortex flow in hy- 
drodynamics, and bacterial chemotaxis in biology [21[|2"2"] ; 
for the present study the theoretical background reviewed 
in Ref. [33J is of particular interest. 

In Ref. [12] the instability of an initially homogeneous 
particle distribution, induced by capillary attraction, has 
been analyzed. The instability is analogous to Jeans' in- 
stability of a self-gravitating fluid 1241125] . except for the 
presence of the screening length A. Since the parameters 
/, A, 7, Rq as well as state variables like particle den- 
sity and temperature are easily tunable, the range of val- 
ues of the parameters for the observation of this insta- 
bility for realistic experimental configurations has been 
explored in Ref. [T2]. Even though the capillary attrac- 
tion between two isolated particles can be relatively weak 
compared with other effects like thermal agitation and 
short-ranged repulsive forces, its long-ranged nature (in 
the sense A Rq) leads to a cumulative effect and to 
an eventual predominance of the capillary force at large 
scales [2B]. Thus, these colloidal monolayers turn out to 
be excellent model systems for the experimental study of 
the collective evolution driven by the pair potential given 
in Eq. |T]). We are aware, however, of only one systematic 
experimental study of this instability with colloidal mono- 
layers |27j , in which the long-ranged capillary interactions 
have been identified as the driving forces of particle ag- 
gregation. 

Here we perform numerical simulations in order to ex- 
tend the results reported in Ref. [H] by going beyond the 
simplifying assumptions underlying the theoretical anal- 
ysis presented there. The paper is organized as follows. 
In Sec. [5] we briefly recall the main theoretical results ob- 
tained in Ref. [12], followed in Sec. [3Jby a description of 
the algorithm of the simulation used here. This includes 
in particular a detailed account of how the long-ranged 
nature of the capillary interaction is handled, similar to 
the approach for simulations of cosmological structure for- 
mation. Section [J presents the analysis of the simulations, 
both for the evolution of an unstable homogeneous distri- 
bution and the evolution of a radially symmetric particle 
density. Particular interest is paid to the dependence of 
the results on the ratio between the capillary length A 
and the size L of the system. This allows one to track 
the crossover between non-integrability (Rq,L <C A) and 



screening (Rq < A « L). Our conclusions are discussed 
and summarized in Sec. [5] 

2 Theoretical model 

In this section we briefly describe the theoretical model 
introduced in Ref. [TJj, which serves as the basis for the 
simulations presented here. 

Upon projection onto the reference plane, for a contact 
angle of = 90°, the configuration of spherical particles 
of radius Rq is mapped onto a configuration of monodis- 
perse discs of radius i?cEL the centers of which form an 
areal number density g(r = (x,y),t), which evolves as a 
function of time under the action of two types of (effective) 
forces: the capillary force and a short-ranged repulsive 
force such as hard core or electrostatic repulsion. Within 
this model, the dynamical evolution of g(r, t) is obtained 
under the following simplifying assumptions. 

• Local equilibrium: The characteristic time scales of 
evolution are sufficiently large so that the forces can be 
described by a local equilibrium ansatz. For the short- 
ranged forces the force per unit area is (V = (d x ,d y )) 

F short (r,t) « -Vp(e= e(r,t)), (2) 

expressed in terms of the equilibrium equation of state for 
the pressure p(g) of the 2D gas of discs in the absence 
of capillary deformations. Isothermal conditions are as- 
sumed. The temperature is a constant parameter fixed by 
the upper and lower fluids acting as thermal baths, and 
thus it is irrelevant for the temporal evolution. 

• Mean field: The capillary force is computed assuming 
that the interfacial deformation is small and adjusts in- 
stantaneously to the momentary particle distribution, i.e., 
it satisfies the Young-Laplace equation for a given particle 
distribution g(r, t). Furthermore, a mean-field approxima- 
tion is used so that the capillary force per unit area can 
be expressed as 

F cap = fgWU, (3) 

with the ensemble averaged mean field interfacial defor- 
mation U(r). By using U (r) one neglects spatial variations 
on small scales, because it is determined self-consistently 
by the mean-field Young-Laplace equation: 

^U-l = - f -g. (4) 

Accordingly it is computed from the average density field 
g(r, t) rather than from the density field conditional to 
the presence of a particle at r [T2] . The pair potential in 
Eq. ([T]) can be expressed as V cav (d) = fUg(d) in terms 
of Green's function Us(d) of Eq. ([4]), i.e., the solution 
for a single-particle source g(r) = S(r). This approxi- 
mation is justified by the long range ~ A (see Eq. ([TJ) 

1 The actual deformation of the interface due to external 
forces / acting on a colloidal particle is small compared with 
its radius Rq [TB]. Therefore, for a contact angle of = 90°, 
the discs have approximately the same radius as the spherical 
colloids. 
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of the capillary force: if g^ is a characteristic homoge- 
neous particle number density, the number of particles 
with which a single particle interacts is effectively of the 
order g h X 2 ~ (A/i? ) 2 > 1. 

• Overdamped motion: Inertia of the particles is ne- 
glected, so that the in-plane velocity field v(r, t) of the 2D 
particle distribution is proportional to the driving force, 



gv(r,t) =r(F sho 



0, 



(5) 



where r, in units of time per mass, is the effective single- 
particle mobility at the interface. For reasons of simplicity 
in Eq. (|SJ) the effect of hydrodynamic interactions is ne- 
glected and r is taken to be constant in space and tirn^E]. 

Under these conditions, mass conservation leads to the 
following form of the continuity equation: 



rV -{Vp- fgVU). 



(6) 



Qualitatively, this equation describes the competition be- 
tween clustering driven by the capillary attraction and 
the homogenization tendency caused by the short-ranged 
forces. The solutions of this equation have been studied 
theoretically in Ref. [H] under certain simplifying condi- 
tions. The corresponding findings can be summarized as 
follows: 

• Jeans' instability: With the homogeneous density 
field one can associate Jeans' time, 



T 



and Jeans' wavenumber, 



K 



7 



/ 2 gjfoh 

7 



(7) 



(8) 



where «h = {ghP'(Qh)) is the isothermal compressibility 
corresponding to the equation of state p(g) for a homoge- 
neous configuration. A linear stability analysis shows that 
a small sinusoidal perturbation gk(r 7 0) = Cgh exp (— zk • r) 
|C| <C 1, of the initially homogeneous distribution g^ 
evolves in time as 

£> k (r, t) = Cg h exp (-ik • r + t/r(k)), (9) 

with a wavenumber-dependent characteristic time 

1 



T(fc) =T 



-1 -1 



[k/KY + (XK)- 



1 



(10) 



From this expression one infers that long-wavelength per- 
turbations with 



k < K c 



1 



(XKf 



(11) 



2 Hydrodynamic interactions play an important role in the 
process of particle trapping by the interface, as studied re- 
cently in Ref. [28]. However, the present theory only considers 
the evolution in the plane of the interface after the particles 
have been trapped there and equilibrium with respect to their 
distribution normal to the interface has been reached. 



grow in time, i.e., r(fc) > 0. Thus the homogeneous state 
is unstable, provided XK > 1. This is known as Jeans' 
instability in a 3D self-gravitating gas }2"41j2"9~] . 

• Cold collapse: This approximation amounts to ne- 
glecting F s hort in Eq. relative to the capillary force. 
(The notion "cold collapse" points to the interpretation of 
this approximation as the zero temperature limit, in which 
the classical pressure p would vanish formally.) This allows 
one to study the nonlinear stage of a capillary-driven col- 
lapse under further simplifying assumptions: considering a 
radially symmetric particle distribution g(r, t) and taking 
A — > oo (which can be phrased as the Newtonian limit), 
one can solve exactly the equation for the radius R{i) of 
a ring of particles. For our present purposes, we consider 
the particular configuration of an initial so-called top-hat 
density profile, 



g(r,0) 



Oh, r < L, 
0, r > L, 



(12) 



with initial size L <C A so that the Newtonian limit holds 
at all times. The trajectory R(t) of any single ring of mat- 
ter can be obtained by specializing the general solution 
given in Ref. |12| to this particular configuration: 



(13) 



expressed in terms of Jeans' time as defined in Eq. ((?]). 
This result predicts that the whole matter collapses at 
the center simultaneously at a time t = T. The ensu- 
ing divergence of this density field is actually preempted 
by the breakdown of the cold-collapse approximation at 
sufficiently high densities, because the force F s h or t halts 
further collapse. 

This approximation has also been used in Ref. |12j 
in order to to predict qualitatively the time of collapse 
of spontaneous fluctuations around a homogeneous back- 
ground density g^: due to the underlying granular nature 
of the density field, in a surface area _4(<C A 2 ) there are un- 
avoidable fluctuations with an amplitude w \J g^A. Within 
the cold-collapse approximation, such initial fluctuations 
would collapse within the time 



7roii ~ Tin yf q^A. 



(14) 



In view of the weak logarithmic dependence, this result 
predicts that for practical purposes fluctuations at all length 
scales collapse on a time scale ~ f. 

In summary, the character of the dynamical evolution 
is set by two length scales: Jeans' length K~ l , above which 
the capillary attraction dominates over short-ranged re- 
pulsion, and the capillary length A, above which the cap- 
illary attraction vanishes. Finally, the short-ranged repul- 
sive forces set the minimum interparticle separation at 
which aggregation stops. 



3 Brownian dynamics simulation 

In this section we describe a theoretical approach for col- 
loidal particles floating at an interface which is distinct 
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from the one given by Eqs. ©- ([!]). To this end we con- 
sider a two-dimensional collection of N disc-like particles 
enclosed in a square box of size Lx L with periodic bound- 
ary conditions. Note that L always denotes the system size 
which is either the side length of a square box or the di- 
ameter of a rotationally symmetric configuration. 

These particles interact via a short-ranged repulsive 
potential V lep (d) which avoids overlap. Specifically, we 
chose a shifted and cut-off Lennard- Jones potential mim- 
icking the hard core repulsion between the particles: 



(or grid) which discretizes the plane. Accordingly, Eq. Q 
leads to 



Vrcpjd) 

knT 



2R 
d 



12 



2Ra 
d 



0, 



d < d r 



d > d r 



(15) 

where R is the radius of the discs and d c /(2R ) = 2 1 / 6 rj 
1.122 the cut-off distance. In addition, a given particle 
is exposed to the capillary force from the other ones (see 
below). As noted before, we neglect any hydrodynamic 
interactions, a systematic treatment of which is beyond 
the scope of the present work. 

Due to the coupling to the two fluid phases, the discs 
experience stochastic forces and their motion is overdamped 
leading to dissipation. This kind of motion is described by 
Brownian dynamics. Applying the Ermak algorithm [30] 
and using the notation introduced in Ref. [31] . the inte- 
gration of the corresponding Langevin equation yields for 
the change of the position of each disc 



rt 



(16) 



where f is a random velocity with (r) = and (f j (i)f j (i')) = 
S(t — t')5ijJ. [31]. The probability distribution of f is taken 
to be uniform with a width w depending on the time step 
At: 

w(At) = V2DAt, D := k B TR (17) 

It turns out that choosing a Gaussian distribution instead 
does not change the results of the present study. The flow 
velocity is v = (r) and the force per particle is f = (F rcp + 

A cap 

)/g (sec Eq. ©). 
This particle-based model constitutes an alternative 
description of the system discussed in the previous sec- 
tion. The random displacements of the colloids in com- 
bination with the action of the short-ranged repulsive 
forces give rise to the macroscopic colloidal pressure p, as 
used in Eq. ([2]), with the corresponding equation of state 
P = p(q) which, in principle, could be determined within 
this model. 

Generically, the capillary force F cap on a single par- 
ticle involves a very large number of neighboring par- 
ticles due to its long range A 3> Ro- Furthermore, we 
arc also interested in the case that A exceeds the side 
length L of the system. Therefore, the capillary inter- 
actions are calculated using the so-called particle-mesh 
algorithm (PM), as applied in cosmological simulations 
which face a similar difficulty (see, e.g., Refs. [32ll33ll34] ). 
The basic idea is to solve Eq. Q via Fourier transfor- 
mation (^[^(r)] = J d 2 r e~ ir ]sL g(r) = g(k)) on a mesh 



1 



so that 



U{r) = T- 



f 



-V-js) U(k) = -ig(k) 



7 



(18a) 



(18b) 



where J- and J 7-1 denote^the Fourier transformation and 
its inverse, respectively. G(k) = — l/(fc 2 + A -2 ) is Green's 
function for this differential equation. The Fourier trans- 
formation can be easily implemented using standard Fast 
Fourier Transform (FFT) routines. The Fourier-transformed 
density held f?(k) is computed from a given g(r) by con- 
structing a discrctized version of the particle density dis- 
tribution. Also G(k) is discretized such that it respects 
the periodic boundary conditions. 

In order to discretize the density, we have employed 
so-called mass assignment schemes using two different ver- 
sions: the so-called Cloud in Cell (CIC) method, in which 
each particle is distributed between four neighboring cells, 
and the so-called Triangular Shaped Cloud (TSC) method 
using nine neighboring cells. The latter one yields an an 
enhanced smoothening of the particle distribution. De- 
pending on the mesh size of the discretization grid, ei- 
ther of these standard methods may be the more accurate 
choice. For a grid with a mesh size (somewhat) smaller 
than the particle size, the distribution to four (CIC) or 
nine (TSC) neighboring cells corresponds to a smeared 
patch the size of which is approximately the actual parti- 
cle size. For the discretized Green's function we use 



G(kh 



1 



- sin 2 (fc x /2) - siii 2 (fc y /2) - A- 2 



with 



_ 2ttZ 

x ~ T' 



2-rtm 



L m G Z. 



(19) 



(20) 



The capillary force is obtained from the interfacial defor- 
mation field via Eq. ([3]): 



>(r) 



f(r) 



V[/C/(r)] 



f 2 - 

t-ikG(k)g(k) 

7 



(21) 



This yields the discrctized capillary force per particle at 
the grid points. In order to obtain the actual force at the 
particle position, one has to interpolate the discretized 
force using the inverse of the chosen mass assignment 
scheme. 

This method offers two important advantages. First, 
it is fast because the number of operations of the PM 
algorithm scales like N g log N g , where A^ is the number 
of nodes of the grid. Second, this method guarantees an 
appropriate incorporation of the long-ranged potential, 
i.e., without introducing any cutoffs together with a high 
accuracy down to small distances of the order of the grid 
spacing r c . We have chosen r c = Rq. 
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In order to save computation time by reducing the 
number of integration steps, and, even more important, 
in order to avoid large displacements Ar by the random 
velocity process, the time step is not fixed during the sim- 
ulation, but it is calculated in each iteration step from 
the requirement that the maximum displacement due to 
the force acting on any individual colloid should not ex- 
ceed i?o/4- This corresponding maximum time step is then 
used in order to calculate the width of the distribution for 
the random displacement according to Eq. ((T7|) . We note 
that this displacement, due to the thermal fluctuations, 
can be larger than the nominal value of i?o/4. This cut- 
off distance appears to be small for simulation boxes with 
L ~ 500i?o- This is, however, a compromise which on one 
hand ensures reasonable computational costs, and on the 
other hand avoids large values of repelling forces. Such 
forces emerge if a particle is displaced randomly on top of 
another particle and encounters the hard core of the re- 
pulsive potential. At such small distances, the r~ 12 term 
in the potential would generate large forces which in turn 
would lead to large velocities, yielding the approximation 
of an overdamped motion and thus the Langevin equation 
for the position inapplicable. 

The structure of the simulation algorithm can be sum- 
marized by the following main steps: 

1. random placement of particles on a plane without over- 
lap; 

2. assignment of particles to a grid by means of a mass 
assignment scheme, resulting in a number density field 
defined on the grid; 

3. Fourier transformation of this density field by applica- 
tion of standard FFT; 

4. computation of the capillary force field on the grid 
according to Eq. (|2"Tj) : 

5. interpolation of this force field to the positions of the 
particles using the inverse mass assignment scheme; 

6. calculation of short-ranged repulsive forces by sum- 
mation over all particles within one box size (minimal 
image convention); 

7. calculation of the maximum possible time step; 

8. integration of the equations of motion by using Er- 
mak's algorithm for Eq. p^]) and loop back to step 
2. 

This algorithm approximates Eq. (O in terms of individual 
particle positions which, when assigned to a density grid, 
represent the approximate evolution of g(r,t). 



4 Results 

In order to probe the predictions of the theoretical model 
reviewed in Sec. [51 for two types of initial configurations 
with periodic boundary conditions we have carried out the 
simulations described above for discs of radius Rq with a 
short-ranged repulsion according to Eq. (|15l) : 



— a rotationally invariant configuration forming a disc of 
radius L surrounded by vacuum, i.e., embedded in an 
empty square box of size L' > 2L. 

We have chosen the temperature to be low enough so 
that these configurations are unstable, i.e., XK > 1 (see 
Eq. HUD ), and the particles therefore tend to lump to- 
gether under the driving capillary force. The effects specif- 
ically due to the long-ranged nature of this interaction 
can be studied by comparing simulation data for various 
values of A. By changing the screening parameter A the 
intcrparticlc attraction can be tuned from being de facto 
short-ranged (A of the order of a few colloid diameters, 
such that the minimum image convention applies) to be- 
ing long-ranged (A — > oo, for which the attraction is for- 
mally analogous to 2D Newtonian gravity). At sufficiently 
low temperatures T, the behavior of the system is deter- 
mined by minimizing its potential energy. Thus, for low T 
the final state for any system of size L is a single cluster 
with maximum close packing of the particles determined 
by the interparticle repulsion (e.g., due to the hard core of 
the discs). This cluster will contain most particles. How- 
ever, the range of the force determines how and on which 
time scale this final state is reached. 

— Rq,L <C A; in this limiting case it is expected that 
phenomena occur which arc similar to the gravitational 
collapse studied in the astrophysical context. In the 
so-called Newtonian limit (A — > oo) the theoretical 
analysis summarized in Sec.Osuggests that Jeans' time 
T (see Eq. ([7])) sets the corresponding characteristic 
time scale of evolution. 

— Ro ^ "C L: in this limit the behavior of the system 
is expected to be similar to that of a common fluid 
and the corresponding mean field description would be 
analogous to a van dcr Waals type equation of state. 
The temporal evolution is expected to follow the sce- 
nario of spinodal decomposition because the initial den- 
sity Qh lies in the unstable region between the liquid 
and vapor coexistence densities. 

One of the goals of the present simulations is to study the 
crossover between these two limiting behaviors. 

We have performed several simulations based on the 
set of parameter values shown in Table [TJ which all cor- 
respond to an experimentally realistic choice of parame- 
ters. Unless stated otherwise, all simulations are based on 
the parameter values quoted in Table [TJ We always cm- 
ploy initial conditions with such diluted densities that the 
equation of state of an ideal gas can be applied in order 
to obtain a rough estimate of Jeans' length from Eqs. (Tj} 
and IE): 



K « —r-= = 2-KQh——. 
7/c B i Kb J 



(22) 



a random homogeneous distribution of non-overlapping 
discs in a square box of size L x L, 



4.1 Random homogeneous distribution as initial 
configuration 

First we address the quantitative prediction in Eqs. © 
and (|10p concerning the early evolution of the initial flue- 
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number of colloids 


N 


2500 


colloid radius 


Ro 


10 /im 


system size 


L 


5 mm 


capillary length 


X 


2.7 mm 


strength of capillary force 


V 


0.89 k B T 


particle mobility 


r 


10.6 x 10~ b s/kg 


mean particle separation 


L/N 1 ' 2 


100 /im 


Jeans' length 


l/K 


43.8 ^m 


Jeans' time 


T 


40780 s 



t = T 



t = 3T 



t = 5T 



Table 1. Default values of the parameters used for all simula- 
tions presented here and certain quantities deduced from them. 
The value of A corresponds to the capillary length of a pure 
water-air interface and the value of F to the effective mobility 
of a sphere half immersed in water at room temperature [12j . 



tuations. In Fig. [5] wc show the growth of the ampli- 
tudes of various modes of the Fourier-transformed ini- 
tial density field, normalized by their values at t = 0. 
The evolution is expected to follow the exponential de- 
pendence given in Eqs. (|9]) and (fT0|) . According to these 
equations, unstable modes are characterized by the con- 
dition k < K c (sec Eq. (fTTjl ) with amplitudes increasing 
with time. We identify the fastest growing mode as the one 
with k = 0.126 K c . Modes with k > K c , the amplitudes of 
which decrease with time, are stable. Up to t ~ 1.5T, the 
agreement with this theoretical prediction is very good. 
Thus the linearized stability analysis is reliable up to the 
characteristic time T of the system. 



* « 



% . 



Fig. 3. Snapshots of the particle distribution in L x L boxes 
with L — 500iio at times T, 3T, and 5T for a "short-ranged" 
capillary attraction with X/L = 0.02, X/L — 0.05, and X/L = 
0.1 (from top to bottom). The number of particles is N = 2500. 
Particles belonging to a cluster are colored in red. 




Fig. 2. Evolution of Fourier modes f?(|k| < K c ,t = 0) of the 
initial density distribution (averaged over shells of constant k — 
|k| and each shifted upwards by 5 for clarity) in a simulation for 
a L x L box, L = 620i?o, with N = 3844 particles. According to 
Eq. (HU) and Table Q] K c = 22.8 mm" 1 and T = 40780s. The 
error bars concerning the statistical errors, as obtained from 
averaging within small time intervals of the order of At « T/40 
and over 20 runs of initial configurations, are smaller than the 
symbol size of the simulation data. The full lines provide the 
corresponding theoretical prediction exp(t/r(k)) as given by 
Eqs. ([9]) and (|10|) . The mean initial number density is = 



The simulations allow us to investigate the highly non- 
linear evolution following the initial clustering stage and 
thus complement the theoretical analysis presented in Ref. [12] 
In Figs. [3] and 2] we show the particle distribution in real 
space at time steps between T and 5T for several values of 
A, while keeping the strength Vq of the capillary potential 
fixed (see Eq. JT])). We define a particle to be a member 
of a cluster (and thus being colored in red) if there are at 
least three neighboring particles at a distance d < 3.25i?o- 
This choice is motivated by the position r m i n ~ 3.25i?o 
of the first minimum in the pair correlation function for a 
dense packing of discs with the repulsion given by Eq. (|15l) , 
as deduced from simulations. For A ~ L (see Fig. 2]) the 
system exhibits the formation of a few clusters on the time 
scale set by Jeans' characteristic time T. These then merge 
quickly leaving the system in its final state of a single clus- 
ter, which contains basically all colloids. For the longest 
range of the capillary potential probed by us, X/L = 10, 
almost all particles in the system collapse into a single 
cluster between t = 3T and t = AT- These observations 
provide a qualitative confirmation of the prediction ac- 
cording to Eq. (fT4"]) that in the Newtonian limit of the 
cold-collapse approximation the time scale of evolution is 
of the order of Jeans' time. 

In order to obtain a more quantitative description of 
these clustering phenomena, the temporal evolution of the 
total number of clusters n c (regardless of their size) is 
presented in Fig. [SJ while Fig. [5] shows the mean cluster 
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t = T 



t = 3T 



t = 5T 



Fig. 4. The same presentation as in Fig.[3]for a "long-ranged" 
capillary attraction with X/L — 0.54, X/L = 3.0, and X/L — 
10.0 (from top to bottom). 



mass defined as 



Efll sn c (s,t) 



SQL) 



(23) 



where n c (s,t) is the number of clusters at time t which 
consist of s particles (note that n c (s < 3, t) = accord- 
ing to the definition of clusters introduced before). As 
A increases, the evolution of these observables converges 
to a A-indcpendcnt limit characterized by scaling time t 
with Jeans' time. This can be identified with the Newto- 
nian limit {X/L — > oo) and with the corresponding onset 
of "long-ranged" dynamics analogous to the gravitational 
collapse known from astrophysical scenarios. In this limit 
the temporal evolution extends over a rather short time 
period (~ 5T), but nevertheless one can interpret Fig. [5] 
qualitatively in terms of a two-stage aggregation process: 
The evolution is initially dominated by the formation of 
new clusters up to a crossover time ~ 1.5T, after which 
the evolution is dominated by the merging of clusters. This 
kind of process is not borne out so obviously by Fig. [6] 
which confirms the steady increase of the average cluster 
mass, and which seems to follow an exponential law be- 
fore saturation, induced by finite-size effects, sets in. This 
latter observation suggests that the average mass S(t) is 
biased towards the largest cluster, so that the early evo- 
lution of the mass si aTge (t) of the largest cluster can be 
described by a simple model of accretion of particles (ei- 
ther isolated or in the form of small clusters) at a rate 
proportional to the capillary attraction by the large clus- 
ter, i.e., 

enlarge 



'o 



CM O 

DC 



'o 



CM O 

rr 




t/T 

Fig. 5. Evolution of the areal number density n c /L 2 of clusters 
for several values of X/L, exhibiting the crossover from "short- 
ranged" attraction (upper plot (a), A < L) to "long-ranged" 
attraction (lower plot (b), A > L). Only representative error 
bars are shown. 




dt 



ks ia 



const. 



(24) 



t/r 

Fig. 6. Evolution of the mean cluster mass S(t) (Eq. 1(230) for 
the same values of A as in Fig. [5] Only representative error 
bars are shown. For X/L < 0.1 the error bars are of the size of 
the symbols. 



As a complement to these observables, we have also 
determined the Minkowski functionals (also known as in- 
trinsic volumes, quermassintcgrale, or curvature integrals) 
associated with the density field. Inter alia, these quanti- 
ties have been applied successfully to the pattern anal- 
ysis of spinodal decomposition in a 2D van-der-Waals 
gas [36] . They provide valuable morphological information 
about the structure formation which is barely accessible 
by other standard probes such as two-point correlation 
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functions because the Minkowski functionals are sensitive 
to higher-order correlations. The properties and applica- 
tions of these quantities in statistical physics are reviewed 
in Ref. [37]. Here we provide only the basic correspond- 
ing informations required for our present purposes. For a 
given geometrical domain T> in the 2D plane, we define 
the three possible Minkowski functionals M„(X>) as (there 
are differences in the literature concerning a conventional 
geometric prefactor) 

M (C) = area of V, (25) 
Mi (V) = length of the boundary of V, (26) 
M 2 (2>) = Euler characteristic of V. (27) 

The Euler characteristic x is a measure of the domain 
connectivity [561157) : 

X = number of clusters minus number of holes, (28) 

where here a "cluster" is defined as a connected subdo- 
main of T> not connected to other "clusters" (i.e., isolated), 
and a "hole" is defined as a "cluster" in the complement 
of T>, i.e., in the domain M. 2 \T). Alternatively, \ can be 
computed also as the integral of the curvature along the 
boundary of T> and thus provides a compact measure of 
whether the boundary of T> is predominantly convex or 
concave. For a given density field g(r) and chosen thresh- 
old values g c , one can construct domains T>(g c ) as the re- 
gions where the density is larger than the threshold value: 



V(Bc) = {r | f?(r) > g c } . 



(29) 



The density field g(r) is computed on a grid with a mass- 
assignment scheme as described in Sec. [3J except that the 
grid is coarser than the one used to compute the cap- 
illary force (for which the grid spacing is r c = R , see 
Table [l}. The purpose of this choice is to enhance the 
smoothening of the density grid; the chosen grid spac- 
ing of ?*Mink. = 5.3i?o offers an acceptable compromise 
between smoothness, accuracy, and computational cost. 
Accordingly, the density field subjected to such a thresh- 
old reduces to a pattern described by pixels and the do- 
main 'D(Qc) is the set of, say, white pixels, for which the 
Minkowski functionals can be computed efficiently (for de- 
tails see Ref. |36j). 

By construction, the Minkowski functionals are func- 
tions of time and of the density threshold g c . The latter 
dependence has been probed in the range between the 
initial number density g^ = (10i?o)~ 2 and the density 
g m = l/(2\/3i?o) ~ 29ph of close packing of discs of ra- 
dius Rq. The qualitative behavior of the functionals docs 
not depend on g c if this value falls within this range (an 
example for the dependence on the threshold value is pre- 
sented in Fig. [TJ. Therefore, in Fig. [S] the time evolution 
of the three Minkowski functionals is shown only for the 
threshold value g c = 0.14p m w Ag^. If the threshold value 
is sufficiently large, as it is the case for g c = Ag^, it is 
expected that in the domain V no "holes" appear, so that 
M?, > and gives directly the number of clusters. Fur- 
thermore, provided the system is ergodic, M is propor- 
tional to the probability distribution of an overdensity, 



i.e., that the density at any point is larger than g c . Mi is 
a direct measure of the length of the "interface" separat- 
ing clusters (i.e., over-dense regions) from the surround- 
ing low-density gas of particles. The quantities M /M 2 , 
Mi/M 2 provide an estimate of the average cluster area 
and perimeter, respectively. We define the "shape factor" 
<fi := Ml / (AttMqM2) which provides a partial, quantita- 
tive characterization of the distribution of cluster shapes 
(see Fig. [5]). It is normalized such that (f> = 1 for a single 
circular cluster. 
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Fig. 7. Evolution of Mi (see Eq. (|26p ) for various values g c 
of the density threshold within the range 0.1 < g c /Qm < 0.8, 
which corresponds to the range 3 < £> c /f?h < 23, for \/L = 
0.54, g m = 29f?h, and Qh = WRo 2 ■ Note that for large times, 
Mi/L converges to the circumference of a close packed disc 
(Mi/L = ^/ N / (tv Qm) w 0.66; not yet visible on the present 
scale) for any density threshold. For g c / g m — 0.1 the error bars 
are twice as large as the symbol, in the other cases they are as 
large as or smaller than the symbol size. 



It is interesting to try to understand the simulation 
results within the framework of spinodal decomposition. 
In our case, the homogeneous initial distribution (which 
is the stable equilibrium state if AA' < 1, see Eqs. ((HJ) 
and (jlip ) is destabilized by a quench into a coexistence 
region between a high-density, liquid phase (g ~ g m ) and 
a rather dilute, gaseous phase. Since the average density 

is substantially smaller than g m , this can be viewed 
as a quench off the center of the coexistence region with 
the liquid state being the minority phase. In this context 
one typically identifies two successive regimes: (i) the early 
stage of phase separation when droplets (i.e., clusters) of 
the minority phase grow out of the initial fluctuation- 
induced seeds by binding individual particles, and (ii) the 
late stage of domain coarsening when the droplets coa- 
lesce. In Ref. , the occurrence of the crossover between 
these two stages has been identified based on the behav- 
ior of the position of the first maximum of the Minkowski 
functionals as function of time. 

Inspection of Figs. HH9] suggests that this spinodal de- 
composition scenario is an appropriate interpretation of 
the present evolution. For X/L = 0.01, which is the small- 
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Fig. 8. Evolution of the Minkowski functionals Mo (a), Mi (b), 
and M 2 (c) (see Eq. (J25J)) for the threshold g c = 4g h = 0.14£> m 
with = WRq 2 and for the same values of A considered in 
Fig. O Mo is proportional to the area of the region where the 
density is higher than the threshold value g c . Mi is propor- 
tional to the contour length of the boundary of this region 
whereas M2 measures the corresponding Euler characteristic. 
All other parameters of the system are chosen as given in Table 
[1] Only representative error bars are shown. 



est value of A studied here, one observes only the initial 
stage of phase separation, marked by a steady increase of 
n c and M2 (number of clusters), S (average mass), Mi 
(accumulated perimeter), and Mq (overdensity area). Ap- 
parently for this small value of A the crossover to domain 
coarsening occurs later than the maximum time of our 
simulations. However, as evidenced for somewhat larger 
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Fig. 9. Same as Fig. [8] for the quantity <j> = Mi/(47rM M 2 ) 
related to the typical shape of clusters and for the threshold 
g c = 0.14g m . Note the different scales in panels (a) and (b). 
Only representative error bars are shown. 



values of A, the shape factor <j> tends to become time- 
independent, indicating a self-similar growth of the clus- 
ters in the domain-coarsening stage. In the opposite limit, 
for the largest values of A explored the crossover to the sec- 
ond stage is clearly observable. It is characterized by the 
occurrence of maxima followed by a steady decrease of the 
number of clusters, of the accumulated perimeter, of the 
overdensity area, and of the shape factor. After the sys- 
tem has collapsed to a single cluster at time t « 6T, the 
increase of the shape factor 0, although subject to fluctu- 
ations, indicates a further compactification and reorgani- 
zation of the clustered particles towards forming a circular 
shape. The average mass S, however, does not exhibit any 
clear feature indicative of the crossover. It eventually sat- 
urates when only a single cluster remains, to which most 
particles are bound. 

We emphasize that this description is valid if the thresh- 
old density g c takes intermediate values. If g c is close to 
the maximum density g m , the signature of the crossover 
is no longer detectable (see Fig. [7]), because in this case 
the Minkowski functionals just probe the very dense re- 
gions, which are formed only at late times and for which 
the dynamics is presumably no longer dominated by cap- 
illary attraction but is instead strongly affected by the 
short-ranged repulsion. 

A characteristic feature of the coarsening stage during 
spinodal decomposition of a common fluid is a self-similar 
dynamical scaling. As we have seen, if A is large the evo- 
lution is too fast and the range of times probed by the 
simulations (a few times 7") is too narrow in order to be 
able to reliably observe a self-similar evolution. However, 
one can still identify the coarsening stage by means of the 
extrema in the Minkowski functionals. This demonstrates 
their robustness even in the extreme case of the almost si- 
multaneous collapse of all length scales, i.e., if most of the 
wavevectors of the system fulfill k < K c (see Eq. (fTTj) ) so 
that the corresponding Fourier components grow rapidly. 
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This supports the applicability of the Minkowski function- 
al as advocated in Ref. [36] . 

In the absence of an elaborated theory, here we can 
provide only a tentative explanation concerning specifi- 
cally the influence of A on the evolution. Our numerical 
observations can be rationalized assuming that the depen- 
dence on A enters only via the two dimensionlcss quanti- 
ties XK and X/L. The data suggest that the Newtonian 
limit (I <C X/L, XK) yields a A-independent behavior 
and that the crossover occurs at a time of the order of 
Jeans' time T, in agreement with Eq. (|14[) . For interme- 
diate values of A (X/L -C 1 <C XK), only the early-time 
evolution is A-independent (in agreement with the ana- 
lytic linear analysis, see Eq. (fTU)) ) whereas the late-time, 
domain-coarsening stage depends significantly on A, but 
the crossover time does so only weakly. Finally, the whole 
evolution depends sensitively on the value of A if it is small 
(X/L <C XK ~ 1): for A > 1/K, the homogeneous state is 
unstable and undergoes spinodal decomposition, whereas 
for A < 1/K the homogeneous state becomes stable as 
estimated from mean-field theory (see Sec. [2]). 

4.2 Radially symmetric initial configuration 

In this subsection we test the cold-collapse approxima- 
tion for radially symmetric configurations. In this context 
"cold" means that the particles move under the action 
of the capillary force only, corresponding to the formal 
limits T — >• (no Brownian motion) and Rq —> (no 
short-ranged repulsion). The first limit is easily attained 
by switching off the random displacement i.e., f = in 
Eq. (|16|) leading to a deterministic Langevin equation (see 
Eq. (|17p). In order to comply with the second limit we 
adopt a sufficiently large initial mean interparticle sepa- 
ration, so that the cold-collapse evolution can be followed 
numerically over some time. This way the breakdown of 
the cold-collapse approximation can be monitored realis- 
tically in those regions where the density approaches its 
largest possible value. All simulations have been carried 
out for N = 451 particles arranged as a disc with initial 
radius L = 183i?o contained in a quadratic simulation box 
with side length L' = 800i?o- Radial symmetry is achieved 
approximately by surrounding the collapsing disc with a 
ringlike stripe of vacuum of thickness w 220i?o- For sim- 
ulations in which the radius of the disc is smaller than A, 
the additional stripe of vacuum ensures that the periodic 
images, which violate the rotational symmetry, contribute 
only weakly to the total force experienced by the particles. 
We have run two sets of simulations: 

1. In the first set, temperature has been switched off 
and particles are initially ordered in suitable concen- 
tric rings with equal spacing (arc length « 12.12i?o) 
along the circumference of all rings, so that each ring 
possesses discrete azimuthal invarianc^l (see inset of 
Fig. [TU|) . The azimuthal position of the particles is pre- 
served during evolution. These somewhat artificial con- 

3 Each ring is azimuthally invariant with respect to multiples 
of the angular separation of nearest neighbors in the same ring. 



ditions provide a way to estimate the influence of the 
finite values of Rq and A on the evolution as compared 
to the theoretical prediction, summarized in Eq. (|13[) . 
2. In the second set the previous conditions are relaxed: 
temperature is switched on as it has been done in the 
simulations discussed in the previous subsection, and 
the initial particle distribution inside the disk is ran- 
dom. With this more realistic study we address the 
effect nonzero temperature has on the theoretical pre- 
diction. 



140 




t/r 

Fig. 10. Cold evolution of the radii R(t) of three rings 
with initial radius Ri/R® = R(t = 0)/Ro = 58, 95, and 135 
for the highly ordered initial configuration of concentric rings 
forming a disc as discussed in the main text (see inset). 
The full lines correspond to Eq. (f30)l with (Aq, Rf/Ro) = 
(1.15, 18), (1.16, 13), and (1.16, 8) top down. The quadratic 
simulation box of size L' = 800i?o contains initially a disc of 
radius L = 183i?o • The number density of particles inside the 
disc is QhRo = 4.28 x 10 -3 The range of the capillary attraction 
is X/L = 1.48 and T = 95340 s. 

For the first set of simulations, Fig. [TU] shows the evo- 
lution of the radial coordinate R(t) for three rings inside 
this disc with different initial radii R{. Motivated by the 
theoretical prediction in Eq. (|13[) , which has been obtained 
in the Newtonian limit within the cold-collapse approxi- 
mation, the time dependence can be fitted rather well by 
the following expression: 



R(t)-R i = (R i -R f )Jl-—. (30) 

This expression differs from Eq. (|13[) in two respects: the 
ring radius R(t) is shifted by subtracting the final radius 
i?f , and the time scale T of collapse is dilated by a factor 
Aq. Deviations at later times, such as the smearing out of 
the square root singularity, can be clearly attributed to the 
non-vanishing size of the particles and the resulting repul- 
sive forces. The effect of the finite compressibility comes 
into play in the late stage of the evolution and prevents the 
vanishing of R(t), so that Rf ^ 0. The factor Aq w 1.16 



J. Bleibel, A. Dommguez, M. Oettel, S. Dietrich: Collective dynamics of colloids at fluid interfaces 



11 



follows from a fit in the range < t < 0.97". However, as 
can be seen in Fig. [TUJ the time AoT offers a reasonable 
estimate for the effective time of collapse in the simula- 
tions 7^oii ~ AqT, defined by the condition R(T co ii) = Rf. 
It turns out that the factor Aq depends only weakly on R{ 
and we conjecture that the delayed effective time of col- 
lapse Tco\h i.e., .Ao > 1, is a consequence of the finite value 
of A. Additionally, residual effects of the tidal forces of the 
periodic images might be non-negligible. Both effects tend 
to weaken the attractive force of the collapsing cluster and 
thus increase T co \\. In order to support this conjecture, in 
Fig. [TT] we show the evolution of the radius R{t) of the 
outermost ring (actually Ri/Ro < L due to binning and 
the applied discretization scheme) for three values of the 
screening length A. Upon increasing A the collapse time 
decreases, eventually leading to a collapse at 7^ n = T, 
i.e., A (X -> oo) — 1 in line with Eq. ([15]). The delay in 
the effective time of collapse T co \i upon decreasing A is 
confirmed in addition by a theoretical calculation of the 
cold-collapse radial trajectories, which is perturbative in 
the parameter 1/A (see Eqs. (4) and (5) in Ref. [26]). This 
effective time must be understood as a fitting parameter 
(according to Eq. ([50)) ) because the analysis presented in 
Ref. [26] indicates that for a small but non vanishing value 
of 1/A the (delayed) collapse does not actually occur at 
the center but at the outer rim of the disk (forming a 
shock wave within the cold-collapse approximation). 
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Fig. 11. Same as Fig.[TU]for the outermost ring (Ri/Ro = 174, 
which is smaller than L due to binning) for X/L — 1, 1.48, 5.9 
and initial disc radius L — 183J?o ■ The full lines correspond to 
the expression in Eq. (|30[) . As expected from Eq. (|13[) it turns 
out that Aq(X — > oo) = 1. 

For the second set of simulations, Fig. Q2] presents the 
results after averaging over the data sets of ten runs cor- 
responding to different realizations of the random initial 
configuration. Equation (|30[) provides a good fit and the 
conclusions following from it hold also for this kind of data 
set. The main differences with Fig. [TUJ are that R(t) ap- 
proaches more smoothly its value Ri of the final, collapsed 
state, and that this happens at later times: the fit of the 



data for the same range < t < 0.9T shows much larger 
deviations at later times so that, although for a fixed value 
of A an approximately common fitting parameter Aq can 
be found, AqT does not provide an equally good estimate 
for the effective time of collapse T co n- These deviations can 
be attributed to the Brownian motion and the random ini- 
tial particle distribution. Both effects tend to spoil the ra- 
dial symmetry and the sole prevalence of capillary attrac- 
tion in driving the evolution and slow down the approach 
to the final state. For the purpose of comparison, the val- 
ues for Jeans' time and wavenumber associated with the 
particle density Qh = N/ (irL 2 ) of the discs are T = 95340 s 
and K = 15.2 jum -1 , respectively. Since A' -1 <C L, L', one 
may have anticipated that the deviations induced by tem- 
perature arc limited to the smallest length scale, i.e., Rf, 
and relevant only after times of order T, as observed in- 
deed. 
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Fig. 12. Same as Fig. 1101 but for an initial random 
distribution of particles (see inset) inside a disc of ini- 
tial radius R(0)/Rq = Ri/Ro and with Brownian motion 
switched on. The radial positions of the particles have been 
averaged over shells with a bin width of Ar = l.hRo- 
The full lines correspond to Eq. (|30j with (Ao, Rf/Ro) = 
(1.16, 18), (1.16, 13), and (1.14, 9) from top to bottom. 



5 Summary and conclusions 

By means of numerical simulations, we have investigated 
the clustering of colloidal particles trapped at a fluid inter- 
face due to capillary attraction. As explained in detail in 
Sec. [21 our numerical approach deals with this long-ranged 
attraction by using a particle-mesh method as employed 
in simulations of self-gravitating fluids. This is motivated 
by the formal analogy between capillary and gravitational 
attraction and because each particle interacts simultane- 
ously with many neighbors. The dynamics is Brownian 
and describes the overdamped motion of particles in flu- 
ids. The key parameters characterizing the macroscopic 
evolution of initial particle configurations are the range 



12 



J. Bleibel, A. Dommguez, M. Oettel, S. Dietrich: Collective dynamics of colloids at fluid interfaces 



A of the capillary attraction (Eq. (JXJ) ) , Jeans' length 1/K 
(Eq. ©), and Jeans' time, T (Eq. ©). 

A distribution of particles is unstable against cluster- 
ing provided the capillary attraction is, in a well-char- 
acterized sense (see Seed]), sufficiently large. Our results 
confirm the predictions derived theoretically in Ref. [T2] : 
the dynamics evolve linearly during the early stage of evo- 
lution up to a time « 1.5T (see Fig. [3]). Also the nonlin- 
ear, radially symmetric collapse in the "cold" (K — > oo) 
as well as in the "Newtonian" (A — > oo) limit is well repro- 
duced (see Figs. [TUHT!?)) . There are small deviations from 
the theoretical predictions which are due to not reaching 
this formal double limit. 

In addition, the simulations provide insights which go 
beyond the simplifying approximations underlying the an- 
alytic theory. In particular, we have addressed the effect 
on the instability of the ratio \/L between the interac- 
tion range A and the system size L. We have analyzed the 
nonlinear evolution by using quantitative characteristics 
which are particularly sensitive to the process of struc- 
ture formation. These are the number and the mass of 
clusters as well as the Minkowski functionals of the den- 
sity field, which have been studied as function of time and 
of X/L (see Figs. [SHE]). In the "short-ranged limit" A<L. 
the evolution proceeds on time scales much larger than T ■ 
During the time probed in our simulations, the evolution 
is dominated by the steady formation of new clusters out 
of the seeds in the fluctuating density field. In the oppo- 
site, "Newtonian" limit A ^> L, a A-indcpcndent behavior 
emerges. The collapse proceeds much faster and is com- 
pleted at a time ~ 5T, thus confirming the conjecture put 
forward in Ref. [12] that in this limit Jeans' time sets the 
temporal scale of the evolution. Despite the speed of this 
evolution, in the clustering process two regimes can be 
identified: an initial stage dominated by the formation of 
new clusters up to a time « 1.57", and a later stage dom- 
inated by cluster merging towards a single, large cluster. 
This is reminiscent of the spinodal decomposition dynam- 
ics in common fluids corresponding to the regimes of phase 
separation and domain coarsening, respectively. 

We interpret this A-dependence of the dynamics as a 
crossover in the behavior from that of a 2D common fluid 
for small A (slow spinodal decomposition) to that of a self- 
gravitating fluid for large A (fast collapse). Our simulation 
results show to which extent the features of each limiting 
case is preserved in the other case and they suggest that 
the dependence on A is smooth. 

Together with the quantitative estimates for T and 
K provided in Ref. Q~2], the present results are expected 
to be helpful for analyzing and interpreting correspond- 
ing experimental data. The simulation algorithm can be 
expanded easily to include the effect of external fields. 
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